Set Up¶
In [1]:
# Import Libraries
import json # Reading and writing JSON
from glob import glob # File pattern matching
import os # Operating system access
import pathlib # File navigation
import holoviews as hv # Interactable plot
import panel as pn # Plot formatting
pn.extension() # Plot formatting
hv.extension('bokeh') # Interactable plot
import earthpy # Work with geospatial data
import xarray as xr # Multi-dimensional arrays
import hvplot.xarray # HVPlot and xarray
import rioxarray as rxr # Raster for xarray
import matplotlib.pyplot as plt # Plotting
import pandas as pd # Work with tabular data
import geopandas as gpd # Wowk with geospatial data
import hvplot.pandas # Plot geospatial data
Download Data¶
The dataset is saved as 'Gila River Vegetation.'
In [2]:
project = earthpy.Project(
'Gila River Vegetation', dirname='gric-vegetation-project')
project.get_data()
The 2020 American Indian Tribal Subdivision National (AITSN) borders from the 2020 US census are included in the data downloaded. More information about the shapefiles can be found here.
In [3]:
# Load in the boundary data
aitsn_gdf = gpd.read_file(project.project_dir / 'tl_2020_us_aitsn')
# Check that it worked
print(type(aitsn_gdf))
aitsn_gdf.head()
<class 'geopandas.geodataframe.GeoDataFrame'>
Out[3]:
| AIANNHCE | TRSUBCE | TRSUBNS | GEOID | NAME | NAMELSAD | LSAD | CLASSFP | MTFCC | FUNCSTAT | ALAND | AWATER | INTPTLAT | INTPTLON | geometry | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 2430 | 653 | 02419073 | 2430653 | Red Valley | Red Valley Chapter | T2 | D7 | G2300 | A | 922036695 | 195247 | +36.6294607 | -109.0550394 | POLYGON ((-109.2827 36.64644, -109.28181 36.65... |
| 1 | 2430 | 665 | 02419077 | 2430665 | Rock Point | Rock Point Chapter | T2 | D7 | G2300 | A | 720360268 | 88806 | +36.6598701 | -109.6166836 | POLYGON ((-109.85922 36.49859, -109.85521 36.5... |
| 2 | 2430 | 675 | 02419081 | 2430675 | Rough Rock | Rough Rock Chapter | T2 | D7 | G2300 | A | 364475668 | 216144 | +36.3976971 | -109.7695183 | POLYGON ((-109.93053 36.40672, -109.92923 36.4... |
| 3 | 2430 | 325 | 02418975 | 2430325 | Indian Wells | Indian Wells Chapter | T2 | D7 | G2300 | A | 717835323 | 133795 | +35.3248534 | -110.0855000 | POLYGON ((-110.24222 35.36327, -110.24215 35.3... |
| 4 | 2430 | 355 | 02418983 | 2430355 | Kayenta | Kayenta Chapter | T2 | D7 | G2300 | A | 1419241065 | 1982848 | +36.6884391 | -110.3045616 | POLYGON ((-110.56817 36.73489, -110.56603 36.7... |
Through some trial and error, I found latitude and longitude boundaries to focus the plot on the Gila River Indian Community. The AIANNHCE value 1310 can then be used to create a new GDF containing just the GRIC Boundaries.
In [4]:
# Generate plot of surrounding area with Tribal locations highlighted
aitsn_gdf.hvplot(
geo=True, tiles='EsriImagery',
frame_width=500,
legend=False, fill_color=None, edge_color='white',
hover_cols='all',
xlim=(-112.4, -111.4),
ylim=(32.8, 33.6)
)
WARNING:param.main: edge_color option not found for polygons plot with bokeh; similar options include: ['muted_color', 'bgcolor', 'line_color']
Out[4]:
In [5]:
# Select and merge the subdivisions we want (1310)
gric_gdf = aitsn_gdf.loc[aitsn_gdf.AIANNHCE=='1310'].dissolve()
# Plot the results with web tile images
gric_gdf.hvplot()
Out[5]:
In [6]:
# Create the boundary GDF
boundary_gdf = aitsn_gdf.loc[aitsn_gdf.AIANNHCE=='1310'].dissolve()
# Plot the results with web tile images
boundary_gdf.hvplot(
geo=True, tiles='EsriImagery',
fill_color='green', line_color='black',
fill_alpha=.3, # must be b/w 0-1
title='Gila River Indian Community',
frame_width=500)
Out[6]:
Set up NDVI data¶
Create and combine a list of the raster files with the unique date for each raster now included as an extra dimension to the (x,y) spatial data.
In [7]:
# Get a sorted list of NDVI tif file paths
ndvi_paths = sorted(list(project.project_dir.rglob('*NDVI*.tif')))
# Display the first and last three files paths to check the pattern
ndvi_paths[:3], ndvi_paths[-3:]
Out[7]:
([PosixPath('/workspaces/data/gric-vegetation-project/gila-river-ndvi/MOD13Q1.061_2001137_to_2022244/MOD13Q1.061__250m_16_days_NDVI_doy2001145000000_aid0001.tif'),
PosixPath('/workspaces/data/gric-vegetation-project/gila-river-ndvi/MOD13Q1.061_2001137_to_2022244/MOD13Q1.061__250m_16_days_NDVI_doy2001161000000_aid0001.tif'),
PosixPath('/workspaces/data/gric-vegetation-project/gila-river-ndvi/MOD13Q1.061_2001137_to_2022244/MOD13Q1.061__250m_16_days_NDVI_doy2001177000000_aid0001.tif')],
[PosixPath('/workspaces/data/gric-vegetation-project/gila-river-ndvi/MOD13Q1.061_2001137_to_2022244/MOD13Q1.061__250m_16_days_NDVI_doy2022209000000_aid0001.tif'),
PosixPath('/workspaces/data/gric-vegetation-project/gila-river-ndvi/MOD13Q1.061_2001137_to_2022244/MOD13Q1.061__250m_16_days_NDVI_doy2022225000000_aid0001.tif'),
PosixPath('/workspaces/data/gric-vegetation-project/gila-river-ndvi/MOD13Q1.061_2001137_to_2022244/MOD13Q1.061__250m_16_days_NDVI_doy2022241000000_aid0001.tif')])
In [8]:
# Loop through each NDVI image to add date dimension
doy_start = -25
doy_end = -19
ndvi_das = []
for ndvi_path in ndvi_paths:
# Get date from file name
doy = ndvi_path.name[doy_start:doy_end]
date = pd.to_datetime(doy, format='%Y%j')
# Open dataset
da = rxr.open_rasterio(ndvi_path, mask_and_scale=True).squeeze()
# Add date dimension and clean up metadata
da = da.assign_coords({'date': date})
da = da.expand_dims({'date': 1})
da.name = 'NDVI'
# Prepare for concatenation
ndvi_das.append(da)
len(ndvi_das)
Out[8]:
154
Combine the list of rasters into a new data xarray "sliced" by the dates.
In [9]:
# Combine NDVI images from all dates
ndvi_da = xr.combine_by_coords(ndvi_das, coords=['date'])
ndvi_da
/tmp/ipykernel_1220/646390900.py:2: FutureWarning: In a future version of xarray the default value for compat will change from compat='no_conflicts' to compat='override'. This is likely to lead to different results when combining overlapping variables with the same name. To opt in to new defaults and get rid of these warnings now use `set_options(use_new_combine_kwarg_defaults=True) or set compat explicitly. ndvi_da = xr.combine_by_coords(ndvi_das, coords=['date']) /tmp/ipykernel_1220/646390900.py:2: FutureWarning: In a future version of xarray the default value for compat will change from compat='no_conflicts' to compat='override'. This is likely to lead to different results when combining overlapping variables with the same name. To opt in to new defaults and get rid of these warnings now use `set_options(use_new_combine_kwarg_defaults=True) or set compat explicitly. ndvi_da = xr.combine_by_coords(ndvi_das, coords=['date'])
Out[9]:
<xarray.Dataset> Size: 48MB
Dimensions: (date: 154, y: 203, x: 382)
Coordinates:
band int64 8B 1
* x (x) float64 3kB -112.3 -112.3 -112.3 ... -111.5 -111.5 -111.5
* y (y) float64 2kB 33.39 33.39 33.38 33.38 ... 32.97 32.97 32.97
spatial_ref int64 8B 0
* date (date) datetime64[ns] 1kB 2001-01-14 2001-01-16 ... 2022-01-24
Data variables:
NDVI (date, y, x) float32 48MB 0.8282 0.6146 ... 0.2146 0.2085In [10]:
# Plot first and last files side by side to visualize progress
first_ndvi = rxr.open_rasterio(ndvi_paths[0], mask_and_scale=True).squeeze()
last_ndvi = rxr.open_rasterio(ndvi_paths[-1], mask_and_scale=True).squeeze()
fig, axes = plt.subplots(1, 2, figsize = (12, 4))
first_ndvi.plot(ax=axes[0], cmap=plt.cm.PiYG, vmin=-1, vmax=1)
axes[0].set_title("Gila River NDVI - 2001")
boundary_gdf.plot(ax=axes[0], edgecolor='black', facecolor='none',
linewidth=1)
last_ndvi.plot(ax=axes[1], cmap=plt.cm.PiYG, vmin=-1, vmax=1)
axes[1].set_title("Gila River NDVI - 2022")
boundary_gdf.plot(ax=axes[1], edgecolor='black', facecolor='none',
linewidth=1)
Out[10]:
<Axes: title={'center': 'Gila River NDVI - 2022'}, xlabel='x', ylabel='y'>
Calculate mean values for NDVI from the beginning of the dataset in 2001 to the water rights settlement in 2004 and an xarray of rasters for each year post 2004.
In [11]:
# Compute the mean NDVI values before water rights were restored
ndvi_pre = (ndvi_da
.sel(date=slice('2001', '2004'))
.mean('date')
.NDVI
)
ndvi_pre.plot()
Out[11]:
<matplotlib.collections.QuadMesh at 0x7c4d6b87b850>
In [12]:
ndvi_post = ndvi_da.sel(date=slice('2005', '2022')).NDVI
ndvi_post_yearly = ndvi_post.groupby('date.year').mean('date')
print (ndvi_post_yearly)
len(ndvi_post_yearly)
<xarray.DataArray 'NDVI' (year: 18, y: 203, x: 382)> Size: 6MB
array([[[0.5380143 , 0.47134283, 0.4491 , ..., 0.17115715,
0.17115715, 0.175 ],
[0.53234285, 0.43379998, 0.44574285, ..., 0.14785714,
0.14785714, 0.16301428],
[0.50931424, 0.4504714 , 0.40725714, ..., 0.1275 ,
0.1056857 , 0.07895714],
...,
[0.18635714, 0.18635714, 0.17112856, ..., 0.17614286,
0.17614286, 0.17889999],
[0.17275715, 0.17275715, 0.16979997, ..., 0.21141426,
0.21141426, 0.2217714 ],
[0.18247142, 0.17312858, 0.18251428, ..., 0.20607142,
0.26507142, 0.2588 ]],
[[0.6633143 , 0.5207143 , 0.4350571 , ..., 0.14278571,
0.14278571, 0.14188571],
[0.6846143 , 0.5030286 , 0.35744286, ..., 0.11975714,
0.11975714, 0.12397142],
[0.5164143 , 0.42664286, 0.3745286 , ..., 0.13028571,
0.09927142, 0.07134286],
...
[0.14237143, 0.14237143, 0.13668571, ..., 0.16885713,
0.16885713, 0.17524286],
[0.1333857 , 0.1333857 , 0.13745713, ..., 0.21079998,
0.21079998, 0.21817143],
[0.15061429, 0.13785715, 0.13837144, ..., 0.16395713,
0.21935715, 0.2031714 ]],
[[0.6115714 , 0.45694283, 0.40044284, ..., 0.16761427,
0.16761427, 0.16799998],
[0.5378714 , 0.44529995, 0.35509998, ..., 0.15104285,
0.15104285, 0.15302856],
[0.33719996, 0.37427142, 0.4122714 , ..., 0.1349 ,
0.11072857, 0.09608571],
...,
[0.12699999, 0.12699999, 0.12062857, ..., 0.15535715,
0.15535715, 0.16347143],
[0.1187857 , 0.1187857 , 0.11717142, ..., 0.19084285,
0.19084285, 0.22217143],
[0.13205715, 0.11748571, 0.12167142, ..., 0.17574285,
0.19661427, 0.21658571]]], shape=(18, 203, 382), dtype=float32)
Coordinates:
band int64 8B 1
* x (x) float64 3kB -112.3 -112.3 -112.3 ... -111.5 -111.5 -111.5
* y (y) float64 2kB 33.39 33.39 33.38 33.38 ... 32.97 32.97 32.97
spatial_ref int64 8B 0
* year (year) int64 144B 2005 2006 2007 2008 ... 2019 2020 2021 2022
Attributes:
units: NDVI
AREA_OR_POINT: Area
Out[12]:
18
In [13]:
# Calculate difference between mean ndvi before water rights and each year following the settlement
ndvi_diff = ndvi_post_yearly - ndvi_pre
In [14]:
# Check dimensions and coordinates for calculated difference
print(ndvi_diff.dims)
print(ndvi_diff.coords)
('year', 'y', 'x')
Coordinates:
band int64 8B 1
* x (x) float64 3kB -112.3 -112.3 -112.3 ... -111.5 -111.5 -111.5
* y (y) float64 2kB 33.39 33.39 33.38 33.38 ... 32.97 32.97 32.97
spatial_ref int64 8B 0
* year (year) int64 144B 2005 2006 2007 2008 ... 2019 2020 2021 2022
In [18]:
# Plot the difference
diff_plot = (
ndvi_diff.hvplot(x='x', y='y', cmap='PiYG', geo=True, groupby='year',
clim=(-0.5, 0.5),
title = 'Gila River Indian Community NDVI\n'
'Comparing 2001-2004 to Following Years')
*
boundary_gdf.hvplot(geo=True, fill_color=None, line_color='black')
)
diff_plot
Out[18]:
BokehModel(combine_events=True, render_bundle={'docs_json': {'9151993b-ad93-44e6-af5e-51ce56d38647': {'version…
Move the slider to the bottom of the plot and set start to 2005.
In [16]:
years = list(diff_plot.kdims[0].values)
years.sort() # just to be safe
diff_plot_bottom = pn.panel(
diff_plot,
widget_location='bottom',
widgets={
'year': pn.widgets.DiscreteSlider(
name='year',
options=years,
value=years[0]
)
}
)
diff_plot_bottom
Out[16]:
BokehModel(combine_events=True, render_bundle={'docs_json': {'2b12ca0f-b559-4be4-afb7-ff4da0790ccc': {'version…
In [19]:
# Save the new plot
diff_plot_bottom.save('ndvi_diff_plot_new.html', embed=True)
WARNING:bokeh.core.validation.check:W-1005 (FIXED_SIZING_MODE): 'fixed' sizing mode requires width and height to be set: figure(id='e0e0ac3a-719d-4a3e-aed1-9ed229c8e554', ...)